This vignette illustrates the use of INLA for spatial prediction using examples from Blangiardo and Cameletti (2015) and Illian, Sørbye, and Rue (2012). For prediction of continuous spatial processes, the Lindgren, Rue, and Lindström (2011) stochastic partial differential equations (SPDE) approach is used to approximate the process through an areal Gaussian Markov random field (GMRF) representation. Finally, Log-Gaussian Cox process models are fit using the pseudodata approach of Simpson et al. (2016).
Blangiardo and Cameletti (2015) section 6.1.
Then the precision matrix \(\mathbf{Q}\) of \(\boldsymbol{\theta}\) is sparse because only neighbors have nonzero coprecisions.
Spatial process \(\xi(\mathbf{s})\).
Finite element approximation…
Toy dataset from Blangiardo and Cameletti (2015).
# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy, pch = 19, asp = 1, main = 'Toy Data')
# Create a mesh for the SPDE method and then plot it.
toy_mesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]), max.edge = c(0.1, 0.2))
plot(toy_mesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
# SPDE projector matrix for estimation.
A_est <- inla.spde.make.A(toy_mesh, as.matrix(SPDEtoy[,c('s1', 's2')]))
# Initialize exponential covariance structure for SPDE.
spde <- inla.spde2.matern(mesh = toy_mesh, alpha = 2)
# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = spde$n.spde)
stack_est <- inla.stack(data = list(y = SPDEtoy$y), A = list(A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')
# Create a grid for prediction.
toy_nx <- 50
toy_ny <- 50
toy_grid <- expand.grid(x = seq(0, 1, length.out = toy_nx), y = seq(0, 1, length.out = toy_ny))
# SPDE projector matrix for prediction.
A_pred <- inla.spde.make.A(mesh = toy_mesh, loc = as.matrix(toy_grid))
# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(y = NA), A = list(A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)
# Fit the model with INLA.
toy_fit <- inla(
y ~ -1 + intercept + f(spatial_field, model = spde),
data = inla.stack.data(stacks),
control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)
# Output posterior summaries.
toy_fit$summary.fixed
toy_fit$summary.hyperpar
# Extract posterior mean of latent spatial field.
index_latent <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- toy_fit$summary.linear.predictor[index_latent, 'mean']
post_sd <- toy_fit$summary.linear.predictor[index_latent, 'sd']
# Plot the posterior mean and SD of the latent spatial field.
plot(im(matrix(post_mean, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior Mean of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
plot(im(matrix(post_sd, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior SD of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005).
# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')
# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
N_QUADS <- 10
QUAD_SIZE <- 50
w_edge <- Frame(bei)$xrange[1]
e_edge <- Frame(bei)$xrange[2]
s_edge <- Frame(bei)$yrange[1]
n_edge <- Frame(bei)$yrange[2]
botleft <- cbind(
runif(N_QUADS, w_edge, e_edge - QUAD_SIZE),
runif(N_QUADS, s_edge, n_edge - QUAD_SIZE)
)
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
cbind(
botleft[r, 1] + c(0, 0, QUAD_SIZE, QUAD_SIZE),
botleft[r, 2] + c(0, QUAD_SIZE, QUAD_SIZE, 0)
)
)})
bei_win <- do.call(
union.owin,
apply(botleft, 1, function(x){return(
owin(x[1] + c(0, QUAD_SIZE), x[2] + c(0, QUAD_SIZE))
)})
)
bei_hole <- bei[complement.owin(bei_win, frame = Frame(bei))]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)
plot(bei_hole, main = 'Observed Region with Holes', pch = '.', cols = 'black')
plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)
# Parameters to experiment with.
MAX_EDGE_LENGTH <- 25
MAX_EDGE_EXT <- 50
MARGIN <- 100
# Mesh covering the site.
bei_boundary <- inla.mesh.segment(loc = do.call(cbind, vertices.owin(Window(bei))))
bei_full_mesh <- inla.mesh.create(
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_LENGTH)
)
plot(bei_full_mesh, asp = 1)
points(bei, pch = 19, cex = 0.25, col = 'red')
# Mesh including a margin outside the site.
margin_mesh <- inla.mesh.2d(
loc = bei_full_mesh$loc[,1:2], # Include nodes from site.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
plot(margin_mesh, asp = 1)
points(bei, pch = 19, cex = 0.25, col = 'red')
# Meshs with coarser resolution in quadrats.
quad_hole <- do.call(
inla.mesh.segment,
lapply(seq_along(bei_interior), function(i){
return(inla.mesh.segment(loc = bei_interior[[i]], grp = i - 1))
})
)
bei_hole_mesh0 <- inla.mesh.create(
boundary = list(bei_boundary, quad_hole),
refine = list(max.edge = MAX_EDGE_LENGTH)
)
plot(bei_hole_mesh0, asp = 1)
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
bei_hole_mesh <- inla.mesh.create(
loc = bei_hole_mesh0$loc[,1:2], # Include nodes from mesh with holes.
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
plot(bei_hole_mesh, asp = 1)
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
# Meshs with finer resolution in quadrats.
quad_bnd <- do.call(
inla.mesh.segment,
lapply(seq_along(bei_interior), function(i){
return(inla.mesh.segment(loc = apply(bei_interior[[i]], 2, rev), grp = i - 1))
})
)
bei_samp_mesh0 <- inla.mesh.create(
boundary = quad_bnd,
refine = list(max.edge = MAX_EDGE_LENGTH)
)
plot(bei_samp_mesh0, asp = 1)
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
bei_samp_mesh <- inla.mesh.create(
loc = bei_samp_mesh0$loc[,1:2], # Include nodes from mesh in quads.
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
plot(bei_samp_mesh, asp = 1)
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
# Meshes with varying resolution in quadrats and a margin.
margin_hole <- inla.mesh.2d(
loc = bei_hole_mesh$loc[,1:2], # Include nodes from mesh with holes.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
plot(margin_hole, asp = 1)
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
margin_samp <- inla.mesh.2d(
loc = bei_samp_mesh$loc[,1:2], # Include nodes from quads.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
plot(margin_samp, asp = 1)
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
NGRID_X <- 40
NGRID_Y <- 20
centers <- gridcenters(
dilation(bei_window_full, max(NGRID_X, NGRID_Y)),
NGRID_X, NGRID_Y)
dx <- sum(unique(centers$x)[1:2] * c(-1, 1)) / 2
dy <- sum(unique(centers$y)[1:2] * c(-1, 1)) / 2
bei_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
message('gridding full: ', system.time(
for(r in seq_len(nrow(bei_df))){
bei_df$count[r] <- sum(bei$x >= bei_df$x[r] - dx &
bei$x < bei_df$x[r] + dx &
bei$y >= bei_df$y[r] - dy &
bei$y < bei_df$y[r] + dy)
bei_df$area[r] <- area(Window(bei)[owin(c(bei_df$x[r] - dx, bei_df$x[r] + dx), c(bei_df$y[r] - dy, bei_df$y[r] + dy))])
}
)['elapsed'], ' seconds')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(bei_df$count, nrow = length(unique(bei_df$x)))), unique(bei_df$x), unique(bei_df$y), unitname = 'meters'), ncolcours = range(bei_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(bei_window_full, border = 'white', add = TRUE)
points(bei, pch = '.', col = 'black')
# SPDE projector matrix for estimation.
full_A_est <- inla.spde.make.A(
bei_full_mesh,
as.matrix(bei_df[bei_df$area > 0, c('x', 'y')])
)
# Initialize exponential covariance structure for SPDE.
full_spde <- inla.spde2.matern(mesh = bei_full_mesh)
# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = full_spde$n.spde)
stack_est <- inla.stack(data = list(count = bei_df$count[bei_df$area > 0], larea = log(bei_df$area[bei_df$area > 0])), A = list(full_A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')
# SPDE projector matrix for prediction.
full_A_pred <- inla.spde.make.A(mesh = bei_full_mesh, loc = as.matrix(bei_df[,c('x', 'y')]))
# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(full_A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(count = NA), A = list(full_A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)
# Fit the model with INLA.
message('gridded full: ', system.time(
bei_full_fit <- inla(
count ~ -1 + intercept + f(spatial_field, model = full_spde),
offset = larea, family = 'poisson',
data = inla.stack.data(stacks),
control.predictor = list(A = inla.stack.A(stacks), compute = TRUE),
verbose = TRUE
)
)['elapsed'], ' seconds')
# Output posterior summaries.
bei_full_fit$summary.fixed
bei_full_fit$summary.hyperpar
# Extract posterior mean of latent spatial field.
index_pred <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- bei_full_fit$summary.linear.predictor[index_pred, 'mean']
post_sd <- bei_full_fit$summary.linear.predictor[index_pred, 'sd']
# Plot the posterior mean and SD of the latent spatial field.
plot(im(t(matrix(post_mean, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior Mean of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')
plot(im(t(matrix(post_sd, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior SD of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')
beihole_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
message('gridding holes: ', system.time(
for(r in seq_len(nrow(beihole_df))){
beihole_df$count[r] <- sum(bei_hole$x >= beihole_df$x[r] - dx &
bei_hole$x < beihole_df$x[r] + dx &
bei_hole$y >= beihole_df$y[r] - dy &
bei_hole$y < beihole_df$y[r] + dy)
beihole_df$area[r] <- area(Window(bei_hole)[owin(c(beihole_df$x[r] - dx, beihole_df$x[r] + dx), c(beihole_df$y[r] - dy, beihole_df$y[r] + dy))])
}
)['elapsed'], ' seconds')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(beihole_df$count, nrow = length(unique(beihole_df$x)))), unique(beihole_df$x), unique(beihole_df$y), unitname = 'meters'), ncolcours = range(beihole_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = '#00000040')
# SPDE projector matrix for estimation.
hole_A_est <- inla.spde.make.A(
bei_hole_mesh0,
as.matrix(bei_df[beihole_df$area > 0, c('x', 'y')])
)
# Initialize exponential covariance structure for SPDE.
hole_spde <- inla.spde2.matern(mesh = bei_hole_mesh0)
beisamp_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
message('gridding sample: ', system.time(
for(r in seq_len(nrow(beisamp_df))){
beisamp_df$count[r] <- sum(bei_samp$x >= beisamp_df$x[r] - dx &
bei_samp$x < beisamp_df$x[r] + dx &
bei_samp$y >= beisamp_df$y[r] - dy &
bei_samp$y < beisamp_df$y[r] + dy)
beisamp_df$area[r] <- area(Window(bei_samp)[owin(c(beisamp_df$x[r] - dx, beisamp_df$x[r] + dx), c(beisamp_df$y[r] - dy, beisamp_df$y[r] + dy))])
}
)['elapsed'], ' seconds')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(beisamp_df$count, nrow = length(unique(beisamp_df$x)))), unique(beisamp_df$x), unique(beisamp_df$y), unitname = 'meters'), ncolcours = range(beisamp_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = '#00000040')
# SPDE projector matrix for estimation.
samp_A_est <- inla.spde.make.A(
bei_samp_mesh0,
as.matrix(bei_df[beisamp_df$area > 0, c('x', 'y')])
)
# Initialize exponential covariance structure for SPDE.
samp_spde <- inla.spde2.matern(mesh = bei_samp_mesh0)
This method relies upon the Lindgren, Rue, and Lindström (2011) approximation of the latent Gaussian field as a linear combination of a finite number of basis functions represented as a GMRF on the nodes of a triangulation of the space. Simpson et al. (2016) use the triangulation for numerical integration of the intensity function and show that the LGCP likelihood factors into the joint distribution of independent Poisson random variables corresponding to the points of the point pattern and the nodes of the triangulation. The model fitting proceeds using INLA to fit a Poisson model to pseudodata.
The pseudodata are constructed as follows.
Then \(y_{i} \sim Poisson(\alpha_{i}\eta_{i})\) where \(\log(\eta_{i})\) is the SPDE representation of the GF at the location of the \(i\)th pseudodatum. See the paper for tedious notation regarding the definition of \(\eta_{i}\). Ultimately, the nodes become Poisson random variables with means equal to the intensity at that their respective locations, observed points become Poisson random variables with means of 1, and the likelihood is approximately proportional to
$ {i=1}^{p+n} {i}^{y_{i}} (-{i} {i}). $
(Is there a missing \(\alpha_{i}\)?)
full_pts <- cbind(bei$x, bei$y)
# Contruct the SPDE A matrix for nodes and points.
full_nV <- bei_full_mesh$n
full_nData <- dim(full_pts)[1]
full_LocationMatrix <- inla.mesh.project(bei_full_mesh, full_pts)$A
full_IntegrationMatrix <- sparseMatrix(i = 1:full_nV, j = 1:full_nV, x = rep(1, full_nV))
full_ObservationMatrix <- rbind(full_IntegrationMatrix, full_LocationMatrix)
# Get the integration weights.
full_IntegrationWeights <- diag(inla.mesh.fem(bei_full_mesh)$c0)
full_E_point_process <- c(full_IntegrationWeights, rep(0, full_nData))
# Create the psuedodata.
full_fake_data <- c(rep(0, full_nV), rep(1, full_nData))
# Fit model to full site.
full_formula <- y ~ -1 + intercept + f(idx, model = full_spde) # No covariates.
full_data <- list(y = full_fake_data, idx = 1:full_nV, intercept = rep(1, full_nV))
message('pseudodata full: ', system.time(
result_full <- inla(
formula = full_formula,
data = full_data,
family = 'poisson',
control.predictor = list(A = full_ObservationMatrix),
E = full_E_point_process,
verbose = TRUE
)
)['elapsed'], ' seconds')
result_full$summary.fixed
result_full$summary.hyperpar
# Fit model to point site with holes.
# Fit model to quadrat-sampled site.
inlabrubei_full_spdf <- as.SpatialPoints.ppp(bei)
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = full_spde) + Intercept
message('inlabru full: ', system.time(
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf, options = list(verbose = TRUE))
)['elapsed'], ' seconds')
bei_full_lgcp$summary.fixed
bei_full_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_full <- predict(bei_full_lgcp, pixels(bei_full_mesh), ~ exp(mySmooth + Intercept))
plot(lambda_full)
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
plot(lambda_full['sd'])
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = hole_spde) + Intercept
message('inlabru with holes: ', system.time(
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf, options = list(verbose = TRUE))
)['elapsed'], ' seconds')
bei_hole_lgcp$summary.fixed
bei_hole_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_hole <- predict(bei_hole_lgcp, pixels(bei_hole_mesh0), ~ exp(mySmooth + Intercept))
plot(lambda_hole)
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
plot(lambda_hole['sd'])
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = samp_spde) + Intercept
message('inlabru quadrats: ', system.time(
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf, options = list(verbose = TRUE))
)['elapsed'], ' seconds')
# Plot posterior means and posterior sd.
lambda_samp <- predict(bei_samp_lgcp, pixels(bei_samp_mesh0), ~ exp(mySmooth + Intercept))
plot(lambda_samp)
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
plot(lambda_samp['sd'])
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.
Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.
Illian, Janine B, Sigrunn H Sørbye, and Håvard Rue. 2012. “A Toolbox for Fitting Complex Spatial Point Process Models Using Integrated Nested Laplace Approximation (Inla).” The Annals of Applied Statistics, 1499–1530.
Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (4): 423–98.
Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.
Simpson, Daniel, Janine B Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue. 2016. “Going Off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.” Biometrika 103 (1): 49–70.